#=
=#


include("../src/qmcmary.jl")
using ..qmcmary

using Random


"""保存hs场的位型"""
function savehs(hscfg, head, fname)
    fil = open(fname, "a")
    write(fil, head)
    write(fil, "\n")
    Nt = size(hscfg)[1]
    for ni in Base.OneTo(Nt)
        write(fil, string(hscfg[ni, :]))
        write(fil, "\n")
    end
    close(fil)
end



function forward(sp, Nt, Nx, Ng, the)
    nx = zeros(Nx)
    ny = sin.(the)
    nz = cos.(the)
    println(ny)
    println(nz)
    fname = "tmp"
    fil = open(fname, "w")
    write(fil, string(the)*"\n")
    close(fil)
    hscfg = Matrix{Int}(undef, Nt, Nx)
    for bi in Base.OneTo(Nt)
        cfg = 2(round.(rand(Nx)) .- 0.5)
        hscfg[bi, :] .= cfg
    end
    sslen, bmats, allbmats, ss = initialize_SS(Nt, Ng, sp, hscfg, nx, ny, nz)
    for _ in Base.OneTo(5)
        hscfg, bmats, allbmats, ss, gf, ph = dqmc_step(Nt, Ng, sp, hscfg, nx, ny, nz, sslen, bmats, allbmats, ss)
    end
    nup_i = zeros(2*Nx)
    tph = 0.0
    nsp = 5
    for idx in Base.OneTo(nsp)
        hscfg, bmats, allbmats, ss, gf, ph = dqmc_step(Nt, Ng, sp, hscfg, nx, ny, nz, sslen, bmats, allbmats, ss)
        tph += ph
        savehs(hscfg, string(idx), fname)
    end
    println("forward sgn: ", tph/nsp)
end


function step_theta1(sp, Nt, Nx, Ng, the)
    fid = open("tmp", "r")
    the2 = readline(fid)
    println(the, the2)
    #
    nx = zeros(Nx)
    ny = sin.(the)
    nz = cos.(the)
    #println(ny)
    #println(nz)
    #
    #重新读取cfg
    ttbar = zeros(Nx)
    tph = 0.0
    bistr = "0"
    tapemap = pbpn_calc_tapemap(sp, nx, ny, nz)
    while !eof(fid)
        bistr = readline(fid)
        #println(bistr)
        hscfg = Matrix{Int}(undef, Nt, Nx)
        for bi in Base.OneTo(Nt)
            cfg = readline(fid)
            cfg = split(cfg[2:end-1], ',')
            for xi in Base.OneTo(Nx)
                hscfg[bi, xi] = strip(cfg[xi]) == "1" ? +1 : -1
            end
        end
        sslen, bmats, allbmats, ss = initialize_SS(Nt, Ng, sp, hscfg, nx, ny, nz)
        ###
        gf, ph = eq_green_scratch(ss)
        #println(ph)
        tph += ph
        #
        nxbar, nybar, nzbar = meas_grad(ss, allbmats, sp, hscfg, nx, ny, nz;
        tapemap=tapemap)
        thetabar = @. nybar * cos(the) -  nzbar * sin(the)
        ttbar += real(thetabar)
    end
    close(fid)
    nsp = parse(Int, bistr)
    println(tph/nsp)
    return ttbar/nsp
end



function optimal_theta1(L)
    #
    the = [0.12166595880333955, -0.12460713674772451, 0.04449671686684168, -0.08196404647845852, 0.10825048664048784, -0.057200788525045114, 0.09953693141895996, -0.07036743487223217]
    println(the)
    Nx = 2*L^2
    hk = lattice_hexagonal(ComplexF64, L, -1.0+0.0im; λ=sqrt(3)+0.0im)
    Ui = 9*ones(Nx)
    Nt = 60
    Ng = 4
    sp = default_splitting(Nt, hk, Ui; Z2=false)
    for ts in Base.OneTo(100)
        Random.seed!(678231)
        forward(sp, Nt, Nx, Ng, the)
        tb = step_theta1(sp, Nt, Nx, Ng, the)
        the += 0.1*tb
    end
end




optimal_theta1(2)

